######## REPLICATION CODE FOR DOST, ENOS, HOCHSCHILD: PSQ 2020, ONLINE APPENDIX A ########

#### SET-UP ####
## set working directory
setwd("/Users/mgdost/Dropbox/race_survey/psq materials/data and replication")
## load packages
library(descr)
library(ggplot2)
library(survey)
library(gridExtra)
library(scales)
library(dotwhisker)
library(psych)
library(WriteXLS)
## load data
# 2017 data set
d17 <- read.csv("cleaned_datasets/harris_dec2017.csv")
# 2019 data set
d19 <- read.csv("cleaned_datasets/harris_feb2019.csv")


#### FIGURE A ####
### 2017
## getting weighted proportions by Trump group
crosstab(d17$QI1_8, d17$stay_leave, prop.c=T, plot=F, weight=d17$weight)
## getting weighted proportions for total sample
freq(d17$QI1_8, w=d17$weight, plot=F)
### 2019
## getting weighted proportions by Trump group
crosstab(d19$QI1R8, d19$stay_leave, prop.c=T, plot=F, weight=d19$weight)
## getting weighted proportions for total sample
freq(d19$QI1R8, w=d19$weight, plot=F)


#### FIGURE B ####
## subset to voters
d19.voters <- subset(d19, stay_leave<=4)
## calculations for table
d1 <- freq(x=d19.voters$male, w=d19.voters$weight, plot=F)[2:1,2]
d2 <- freq(x=d19.voters$racerec, w=d19.voters$weight, plot=F)[1:4,2]
d3 <- freq(x=d19.voters$agerec, w=d19.voters$weight, plot=F)[1:4,2]
d4 <- freq(x=d19.voters$educ, w=d19.voters$weight, plot=F)[1:4,2]
d5 <- freq(x=d19.voters$increc_4, w=d19.voters$weight, plot=F)[1:4,2]
d6 <- freq(x=d19.voters$ideo, w=d19.voters$weight, plot=F)[1:5,2]
d7 <- freq(x=d19.voters$ptystr, w=d19.voters$weight, plot=F)[1:7,2]
d8 <- freq(x=d19.voters$employ, w=d19.voters$weight, plot=F)[1:4,2]
d9 <- freq(x=d19.voters$QD16, w=d19.voters$weight, plot=F)[1:3,2]
d10 <- freq(x=d19$howvote, w=d19$weight, plot=F)[1:4,3]
# combining results and adding labels
d <- cbind.data.frame(c(d1,d2,d3,d4,d5,d6,d7,d8,d9,d10),
                      c("Men","Women","White NH","Black NH","Hispanic","Other NH",
                        "18-29","30-49","50-64","65+","HS or less","Some college","College grad","Post grad",
                        "Income <$35K","Income $35-75K","Income $75-125K","Income $125K+",
                        "Strong Conserv","Lean Conserv","Moderate","Lean Liberal","Strong Liberal",
                        "Strong Rep","Not Strong Rep","Lean Rep","Pure Ind","Lean Dem","Not Strong Dem","Strong Dem",
                        "Employed full-time","Employed part-time","Not employed, not looking","Not employed, looking",
                        "Urban","Suburban","Rural","Voted Trump","Clinton","Other","Voted; ref vote choice"))
names(d) <- c("freq","demo")
d.rot <- cbind.data.frame(d$demo,round(d$freq))
write.table(d.rot, file = "online_appendix_a/plots/FigB_demo_table.txt", sep = ",", quote = FALSE, row.names = F)
rm(d,d.rot,d19.voters,d1,d2,d3,d4,d5,d6,d7,d8,d9,d10)


#### FIGURE C ####
### part 1: cosmopolitanism
d19$COSMOPOLITANISMR4_rec <- recode(d19$COSMOPOLITANISMR4, "1=5;2=4;3=3;4=2;5=1")
d19$COSMOPOLITANISMR6_rec <- recode(d19$COSMOPOLITANISMR6, "1=5;2=4;3=3;4=2;5=1")
cosmo.vars <- c("COSMOPOLITANISMR1", "COSMOPOLITANISMR2","COSMOPOLITANISMR3", "COSMOPOLITANISMR4_rec", "COSMOPOLITANISMR5", "COSMOPOLITANISMR6_rec")
cosmo.cor <- data.frame(cor(d19[,cosmo.vars]))
eigen(cor(d19[,cosmo.vars]))$values # keep two components
cosmo.comp <- principal(d19[,cosmo.vars], nfactors = 2, rotate= "varimax")
cosmo.comp.loadings <- rbind.data.frame(unclass(cosmo.comp$loadings), cosmo.comp$Vaccounted)
### part 2: is it racist?
racism.vars <- c("RACISM_OKR1", "RACISM_OKR2", "RACISM_OKR3", "RACISM_OKR4", "RACISM_OKR5","RACISM_OKR6", "RACISM_OKR7", "RACISM_OKR8", "RACISM_OKR9", "RACISM_OKR10")
racism.cor <- data.frame(cor(d19[,racism.vars]))
eigen(cor(d19[,racism.vars]))$values # keep two components
racism.comp <- principal(d19[,racism.vars], nfactors = 2, rotate= "varimax")
racism.comp.loadings <- rbind.data.frame(unclass(racism.comp$loadings), racism.comp$Vaccounted)
### part 3: truly American
american.vars <- c("TRULY_AMERICANR1", "TRULY_AMERICANR2", "TRULY_AMERICANR3", "TRULY_AMERICANR4", "TRULY_AMERICANR5", "TRULY_AMERICANR6", "TRULY_AMERICANR7", "TRULY_AMERICANR8")
american.cor <- data.frame(cor(d19[,american.vars]))
eigen(cor(d19[,american.vars]))$values # keep two components
american.comp <- principal(d19[,american.vars], nfactors = 2, rotate= "varimax")
american.comp.loadings <- rbind.data.frame(unclass(american.comp$loadings), american.comp$Vaccounted)
### part 4: law favorability
legis.vars <- c("LEGISLATIONR1","LEGISLATIONR2","LEGISLATIONR3","LEGISLATIONR4","LEGISLATIONR5","LEGISLATIONR6","LEGISLATIONR7","LEGISLATIONR8","LEGISLATIONR9","LEGISLATIONR10","LEGISLATIONR11","LEGISLATIONR12","LEGISLATIONR13", "LEGISLATIONR14", "LEGISLATIONR15")
legis.cor <- data.frame(cor(d19[,legis.vars]))
eigen(cor(d19[,legis.vars]))$values # keep three components
legis.comp <- principal(d19[,legis.vars], nfactors = 3, rotate= "varimax")
legis.comp.loadings <- rbind.data.frame(unclass(legis.comp$loadings), legis.comp$Vaccounted)
### part 5: tax favorability
taxes.vars <- c("TAXESR1","TAXESR2","TAXESR3","TAXESR4","TAXESR5","TAXESR6","TAXESR7","TAXESR8","TAXESR9","TAXESR10","TAXESR11","TAXESR12","TAXESR13", "TAXESR14", "TAXESR15")
taxes.cor <- data.frame(cor(d19[,taxes.vars]))
eigen(cor(d19[,taxes.vars]))$values # keep four components
taxes.comp <- principal(d19[,taxes.vars], nfactors = 4, rotate= "varimax")
taxes.comp.loadings <- rbind.data.frame(unclass(taxes.comp$loadings), taxes.comp$Vaccounted)
### saving out correlation matrices and principal component analyses
WriteXLS(list(cosmo.cor,cosmo.comp.loadings,racism.cor,racism.comp.loadings,american.cor,american.comp.loadings,
              legis.cor,legis.comp.loadings,taxes.cor,taxes.comp.loadings),"online_appendix_a/FigC_scale_analyses_2019.xlsx",
         SheetNames=c("cosmo_CORR","cosmo_PCA","racism_ok_CORR", "racism_ok_PCA","american_CORR","american_PCA",
                      "legislation_CORR","legislation_PCA","taxes_CORR","taxes_PCA"),row.names = T)


#### FIGURE D ####
## calculations prior to figure creation 
t <- crosstab(d19$ideo, d19$stay_leave, prop.c=T, plot=F, weight=d19$weight)
t2 <- crosstab(d19$ptystr, d19$stay_leave, prop.c=T, plot=F, weight=d19$weight)
t3 <- crosstab(d19$agerec, d19$stay_leave, prop.c=T, plot=F, weight=d19$weight)
t4 <- crosstab(d19$increc_4, d19$stay_leave, prop.c=T, plot=F, weight=d19$weight)
t5 <- crosstab(d19$male, d19$stay_leave, prop.c=T, plot=F, weight=d19$weight)
t6 <- crosstab(d19$racerec, d19$stay_leave, prop.c=T, plot=F, weight=d19$weight)
t7 <- crosstab(d19$educ, d19$stay_leave, prop.c=T, plot=F, weight=d19$weight)
t8 <- crosstab(d19$employ, d19$stay_leave, prop.c=T, plot=F, weight=d19$weight)
t9 <- crosstab(d19$QD16, d19$stay_leave, prop.c=T, plot=F, weight=d19$weight) # with QD16 values: 1=urban;2=suburban;3=rural
df1 <- data.frame((t$prop.col)*100)
names(df1) <- c("demo","TG","TG_freq")
df1$demo <- car::recode(df1$demo, "1='Strong Liberal';2='Lean Liberal';3='Moderate';4='Lean Conserv';5='Strong Conserv'")
df1$col <- 8
df2 <- data.frame((t2$prop.col)*100)
names(df2) <- c("demo","TG","TG_freq")
df2$demo <- car::recode(df2$demo, "1='Strong Rep';2='Not Strong Rep';3='Lean Rep';4='Pure Ind';5='Lean Dem';6='Not Strong Dem';7='Strong Dem'")
df2$col <- 9
df3 <- data.frame((t3$prop.col)*100)
names(df3) <- c("demo","TG","TG_freq")
df3$demo <- car::recode(df3$demo, "1='18-29';2='30-49';3='50-64';4='65+'")
df3$col <- 3
df4 <- data.frame((t4$prop.col)*100)
names(df4) <- c("demo","TG","TG_freq")
df4$demo <- car::recode(df4$demo, "1='Income <$35K';2='Income $35-75K';3='Income $75-125K';4='Income $125K+'")
df4$col <- 5
df5 <- data.frame((t5$prop.col)*100)
names(df5) <- c("demo","TG","TG_freq")
df5$demo <- car::recode(df5$demo, "1='Men';0='Women'")
df5$col <- 1
df6 <- data.frame((t6$prop.col)*100)
names(df6) <- c("demo","TG","TG_freq")
df6$demo <- car::recode(df6$demo, "1='White NH';2='Black NH';3='Hispanic';4='Other NH'")
df6$col <- 2
df7 <- data.frame((t7$prop.col)*100)
names(df7) <- c("demo","TG","TG_freq")
df7$demo <- car::recode(df7$demo, "1='HS or less';2='Some college';3='College grad';4='Post grad'")
df7$col <- 4
df8 <- data.frame((t8$prop.col)*100)
names(df8) <- c("demo","TG","TG_freq")
df8$demo <- car::recode(df8$demo, "1='Employed full-time';2='Employed part-time';3='Not employed, not looking';4='Not employed, looking'")
df8$col <- 6
df9 <- data.frame((t9$prop.col)*100)
names(df9) <- c("demo","TG","TG_freq")
df9$demo <- car::recode(df9$demo, "1='Urban';2='Suburban';3='Rural'")
df9$col <- 7
df <- rbind.data.frame(df1,df2,df3,df4,df5,df6,df7,df8,df9)
df$TG <- as.numeric(as.character(df$TG))
df <- subset(df, TG<=4)
df$demo <- factor(df$demo,levels = c("Men","Women",paste(unique(df6$demo)),paste(unique(df3$demo)),
                                     paste(unique(df7$demo)),paste(unique(df4$demo)),
                                     paste(unique(df8$demo)),paste(unique(df9$demo)),"Strong Conserv",
                                     "Lean Conserv","Moderate","Lean Liberal","Strong Liberal",
                                     paste(unique(df2$demo))))
rm(df1,df2,df3,df4,df5,df6,df7,df8,df9,t,t2,t3,t4,t5,t6,t7,t8,t9)
# getting standard errors
dsgn <- svydesign(ids = ~RECORD, weights = ~weight, data = subset(d19, stay_leave<=4)) 
df.se <- data.frame(matrix(c(rep(1,4),rep(2,4),rep(3,4),rep(4,4),rep(5,4),rep(6,4),rep(7,4),rep(8,4),rep(9,4),
                             rep(1:4,9),rep(NA,36)),ncol=3,nrow=length(unique(df$col))*4))
names(df.se) <- c("col","TG","sd")
df.se[1:4,3] <- svyby(~male, ~stay_leave,dsgn, svymean)[,3]
df.se[5:8,3] <- svyby(~racerec, ~stay_leave,dsgn, svymean)[,3]
df.se[9:12,3] <- svyby(~agerec, ~stay_leave,dsgn, svymean)[,3]
df.se[13:16,3] <- svyby(~educ, ~stay_leave,dsgn, svymean)[,3]
df.se[17:20,3] <- svyby(~increc_4, ~stay_leave,dsgn, svymean,na.rm=T)[,3]
df.se[21:24,3] <- svyby(~employ, ~stay_leave,dsgn, svymean,na.rm=T)[,3]
df.se[25:28,3] <- svyby(~QD16, ~stay_leave,dsgn, svymean,na.rm=T)[,3]
df.se[29:32,3] <- svyby(~ideo, ~stay_leave,dsgn, svymean)[,3]
df.se[33:36,3] <- svyby(~ptystr, ~stay_leave,dsgn, svymean,na.rm=T)[,3]
df <- merge(df, df.se, by = c("col","TG"))
df$sd <- df$sd*100
# subsetting into groups
df_stay <- subset(df, TG==1)
df_leave <- subset(df, TG==2)
df_join <- subset(df, TG==3)
df_never <- subset(df, TG==4)
rm(df, df.se)
### creating individual plots for each Loyalist/Switcher group for Figure C
my.palette <- c("#8DD3C7","#FFFFB3","#BEBADA","#FB8072","#80B1D3","#FDB462","#B3DE69","#FF99CC","#BC80BD")
## Republican Loyalists plot
p1 <- ggplot(df_stay, aes(x=factor(demo, levels = rev(levels(factor(demo)))), y=TG_freq, fill = factor(col))) + 
  geom_bar(stat = "identity") + geom_text(aes(y=TG_freq/2, label=round(TG_freq)), color="black", size=2.5) +
  geom_errorbar(aes(ymin=TG_freq-sd, ymax=TG_freq+sd), width=.5) +
  scale_fill_manual(values = my.palette) + theme(legend.position="none",axis.ticks.y = element_blank()) + ylab("tk") + xlab("") +
  ggtitle("Republican Loyalists") + coord_flip() + theme(plot.title = element_text(size=10), axis.text = element_text(size=8), axis.title = element_text(color="white")) +
  scale_x_discrete(limits = c(levels(df_stay$demo)[37:31],"b8",levels(df_stay$demo)[30:26],"b7",
                              levels(df_stay$demo)[25:23],"b6",levels(df_stay$demo)[22:19],"b5",
                              levels(df_stay$demo)[18:15],"b4",levels(df_stay$demo)[14:11],"b3",
                              levels(df_stay$demo)[10:7],"b2",levels(df_stay$demo)[6:3],"b1",
                              levels(df_stay$demo)[2:1]),
                   labels=c("b8"="","b7"="","b6"="","b5"="","b4"="","b3"="","b2"="","b1"="")) +
  scale_y_continuous(breaks=c(0,25,50,75), limits = c(-23,90)) 
## Republican Switchers plot
p2 <- ggplot(df_leave, aes(x=factor(demo, levels = rev(levels(factor(demo)))), y=TG_freq, fill = factor(col))) + 
  geom_bar(stat = "identity") + geom_text(aes(y=TG_freq/2, label=round(TG_freq)), color="black", size=2.5) +
  geom_errorbar(aes(ymin=TG_freq-sd, ymax=TG_freq+sd), width=.5) +
  scale_fill_manual(values = my.palette) + theme(legend.position="none",axis.ticks.y = element_blank()) + ylab("tk") + xlab("") +
  ggtitle("Republican Switchers") + coord_flip() + theme(plot.title = element_text(size=10), axis.text = element_text(size=8), axis.title = element_text(color="white")) +
  scale_x_discrete(limits = c(levels(df_stay$demo)[37:31],"b8",levels(df_stay$demo)[30:26],"b7",
                              levels(df_stay$demo)[25:23],"b6",levels(df_stay$demo)[22:19],"b5",
                              levels(df_stay$demo)[18:15],"b4",levels(df_stay$demo)[14:11],"b3",
                              levels(df_stay$demo)[10:7],"b2",levels(df_stay$demo)[6:3],"b1",
                              levels(df_stay$demo)[2:1]),
                   labels=c("b8"="","b7"="","b6"="","b5"="","b4"="","b3"="","b2"="","b1"="")) +
  scale_y_continuous(breaks=c(0,25,50,75), limits = c(-23,90)) 
## Democratic Switchers plot
p3 <- ggplot(df_join, aes(x=factor(demo, levels = rev(levels(factor(demo)))), y=TG_freq, fill = factor(col))) + 
  geom_bar(stat = "identity") + geom_text(aes(y=TG_freq/2, label=round(TG_freq)), color="black", size=2.5) +
  geom_errorbar(aes(ymin=TG_freq-sd, ymax=TG_freq+sd), width=.5) +
  scale_fill_manual(values = my.palette) + theme(legend.position="none",axis.ticks.y = element_blank()) + ylab("% identifying as each") + xlab("") +
  ggtitle("Democratic Switchers") + coord_flip() + theme(plot.title = element_text(size=10), axis.text = element_text(size=8)) +
  scale_x_discrete(limits = c(levels(df_stay$demo)[37:31],"b8",levels(df_stay$demo)[30:26],"b7",
                              levels(df_stay$demo)[25:23],"b6",levels(df_stay$demo)[22:19],"b5",
                              levels(df_stay$demo)[18:15],"b4",levels(df_stay$demo)[14:11],"b3",
                              levels(df_stay$demo)[10:7],"b2",levels(df_stay$demo)[6:3],"b1",
                              levels(df_stay$demo)[2:1]),
                   labels=c("b8"="","b7"="","b6"="","b5"="","b4"="","b3"="","b2"="","b1"="")) +
  scale_y_continuous(breaks=c(0,25,50,75), limits = c(-23,90)) 
## Democratic Loyalists plot
p4 <- ggplot(df_never, aes(x=factor(demo, levels = rev(levels(factor(demo)))), y=TG_freq, fill = factor(col))) + 
  geom_bar(stat = "identity") + geom_text(aes(y=TG_freq/2, label=round(TG_freq)), color="black", size=2.5) +
  geom_errorbar(aes(ymin=TG_freq-sd, ymax=TG_freq+sd), width=.5) +
  scale_fill_manual(values = my.palette) + theme(legend.position="none",axis.ticks.y = element_blank()) + ylab("% identifying as each") + xlab("") +
  ggtitle("Democratic Loyalists") + coord_flip() + theme(plot.title = element_text(size=10), axis.text = element_text(size=8)) +
  scale_x_discrete(limits = c(levels(df_stay$demo)[37:31],"b8",levels(df_stay$demo)[30:26],"b7",
                              levels(df_stay$demo)[25:23],"b6",levels(df_stay$demo)[22:19],"b5",
                              levels(df_stay$demo)[18:15],"b4",levels(df_stay$demo)[14:11],"b3",
                              levels(df_stay$demo)[10:7],"b2",levels(df_stay$demo)[6:3],"b1",
                              levels(df_stay$demo)[2:1]),
                   labels=c("b8"="","b7"="","b6"="","b5"="","b4"="","b3"="","b2"="","b1"="")) +
  scale_y_continuous(breaks=c(0,25,50,75), limits = c(-23,90)) 
## breakdown of our sample 
d19sub <- subset(d19, stay_leave<=4)
t0 <- descr::freq(x=d19sub$stay_leave, w=d19sub$weight)[1:4,2]
df0 <- data.frame(matrix(c(4,2,3,1,t0,0,0,0,0),nrow=4,ncol=3))
names(df0) <- c("TG","TG_freq","x")
df0$TG <- as.factor(df0$TG)
df0$label_loc <- NA
df0$label_loc[df0$TG==1] <- df0$TG_freq[df0$TG==1]/2
df0$label_loc[df0$TG==3] <- df0$TG_freq[df0$TG==3]/2 + df0$TG_freq[df0$TG==1]
df0$label_loc[df0$TG==2] <- df0$TG_freq[df0$TG==2]/2 + df0$TG_freq[df0$TG==1] + df0$TG_freq[df0$TG==3]
df0$label_loc[df0$TG==4] <- df0$TG_freq[df0$TG==4]/2 + df0$TG_freq[df0$TG==1] + df0$TG_freq[df0$TG==2] + df0$TG_freq[df0$TG==3]
df0$N <- as.integer(table(d19sub$stay_leave))
df0$groupname <- c("Rep Loyalists","Rep\nSwitch","Dem\nSwitch","Dem Loyalists")
p5 <- ggplot(data=df0, aes(x=x, y=TG_freq, fill=factor(TG, levels = c("4","2","3","1")))) +
  geom_bar(stat="identity") +
  geom_text(aes(x=0.275,y=label_loc, label=groupname), color="black", size=2, fontface="bold") +
  geom_text(aes(x=-0.1,y=label_loc, label=paste(round(TG_freq),"%",sep="")), color="white", size=2) +
  geom_text(aes(x=-0.35,y=label_loc, label=N), color="white", size=2) +
  geom_text(aes(x=-0.35,y=-1, label="Unwtd N"), hjust=1,color="black", size=2) +
  geom_text(aes(x=-0.1,y=-1, label="Proportion"), hjust=1,color="black", size=2) +
  xlab("") + ylab("") + ylim(-9,100) +
  theme(legend.position="none",plot.title = element_text(hjust=0.5),
        panel.background = element_rect(fill = "white",colour = "white"),
        panel.border = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text.x = element_blank(), axis.ticks = element_blank(),
        axis.text.y = element_blank()) + 
  coord_flip() + scale_fill_manual(values = c("darkorange", alpha("darkorange", .45), alpha("purple3", .45), "purple3"))
### creating final plot of Figure C
p <- grid.arrange(p5, p1, p2, p4, p3, ncol=2,
                  layout_matrix = rbind(c(1,1),c(2,3),c(4,5)),
                  widths = c(2.5,2.5), heights = c(.4,2.6,2.75))
ggsave("online_appendix_a/plots/FigD_GroupProfile2019_CIs.png", p, width = 8, height = 10.5, units="in")
rm(d19sub,df_never,df_join,df_leave,df_stay,df0,dsgn,my.palette,p,p1,p2,p3,p4,p5,t0)


#### FIGURE E ####
## see separate Rmd file


#### PREP FOR FIGURE F, parts 1a-5b ####
## setting up datasets for plots
# subsetting datasets and standardizing scales
vars <- c("cosmo_scale2","american_scale_v1","american_scale_v2","taxes_scale_v1","laws_toofavorable","laws_toofavorable_wealthycorps", "trump_immig_app", #nativism
          "ECON_LADDER_past_future","ECON_LADDER_past_present","ECON_LADDER_present_future","personalfinance","PARENTS_rec","CHILDREN_rec", #economic anxiety (egocentric)
          "econ_wrongtrack","econ_strength", "country_wrongtrack","trump_disapp_econ","trump_disapp_stimjobs", #economic anxiety (sociotropic)
          "BREAK_RULES_rec","US_DARK_rec","COMPROMISE","taxes_scale_v2", #authoritarisnism
          "racial_resent_scale", # racism
          "TAXESR1","TAXESR2","TAXESR3","TAXESR4","TAXESR5","TAXESR6","TAXESR7","TAXESR8","TAXESR9","TAXESR10",
          "TAXESR11","TAXESR12","TAXESR13","TAXESR14","TAXESR15","LEGISLATIONR1","LEGISLATIONR2","LEGISLATIONR3",
          "LEGISLATIONR4","LEGISLATIONR5","LEGISLATIONR6","LEGISLATIONR7","LEGISLATIONR8","LEGISLATIONR9",             
          "LEGISLATIONR10","LEGISLATIONR11","LEGISLATIONR12","LEGISLATIONR13","LEGISLATIONR14","LEGISLATIONR15",
          "stay_leave","weight","RECORD",
          "ptystr")
d19.sub <- d19[vars]
x <- data.frame(sapply(c(1:53), function(i){scales::rescale(as.numeric(unlist(d19.sub[,i])), to = c(-1,1))}))
d19.sub <- cbind.data.frame(x,d19.sub[,54:57])
colnames(d19.sub) <- vars
# subsetting among partisans
d19rep <- subset(d19.sub, ptystr<4)
d19dem <- subset(d19.sub, ptystr>4)

#### RUN THE FOLLOWING CODE FOR DEMOCRAT SUBSET FIRST, THEN REPUBLICAN #### 
#### TO REPLICATE FIGURE E, parts 1a-5b
### this way, the plotted items appear in the same order for Dem and Rep plots
### run only one of the next two lines!
d19.sub <- d19dem
#d19.sub <- d19rep

# subset datasets
d19_St_Lv <- subset(d19.sub, stay_leave==1 | stay_leave==2)
d19_Jo_Nv <- subset(d19.sub, stay_leave==3 | stay_leave==4)
# setting up survey designs
dsgn.stay19 <- svydesign(ids = ~RECORD, weights = ~weight, data = d19_St_Lv)
dsgn.never19 <- svydesign(ids = ~RECORD, weights = ~weight, data = d19_Jo_Nv)
# prep for plot DFs
x <- c("Rep Loyalists","Rep Switchers","Dem Switchers","Dem Loyalists")

#### FIGURE F - part 1 ####
## doing calculations for plots
# econ_wrongtrack
dsgn.stay19.wr <- subset(dsgn.stay19, !is.na(dsgn.stay19$variables$econ_wrongtrack)) # dealing w/NAs
dsgn.never19.wr <- subset(dsgn.never19, !is.na(dsgn.never19$variables$econ_wrongtrack)) # dealing w/NAs
y2 <- c(svyby(~econ_wrongtrack, ~ stay_leave, dsgn.stay19.wr, svymean)$econ_wrongtrack,svyby(~econ_wrongtrack, ~ stay_leave, dsgn.never19.wr, svymean)$econ_wrongtrack)
y2_error <- c(sqrt(svyby(~econ_wrongtrack, ~stay_leave, dsgn.stay19.wr, svyvar)$econ_wrongtrack)/sqrt(nrow(dsgn.stay19.wr)),sqrt(svyby(~econ_wrongtrack, ~stay_leave, dsgn.never19.wr, svyvar)$econ_wrongtrack)/sqrt(nrow(dsgn.never19.wr)))
# econ_strength 
y3 <- c(svyby(~econ_strength, ~ stay_leave, dsgn.stay19, svymean)$econ_strength,svyby(~econ_strength, ~ stay_leave, dsgn.never19, svymean)$econ_strength)
y3_error <- c(sqrt(svyby(~econ_strength, ~stay_leave, dsgn.stay19, svyvar)$econ_strength)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~econ_strength, ~stay_leave, dsgn.never19, svyvar)$econ_strength)/sqrt(nrow(d19_Jo_Nv)))
y3 <- y3*-1 # reversing scale
# country_wrongtrack 
dsgn.stay19.tr <- subset(dsgn.stay19, !is.na(dsgn.stay19$variables$country_wrongtrack)) # dealing w/NAs
dsgn.never19.tr <- subset(dsgn.never19, !is.na(dsgn.never19$variables$country_wrongtrack)) # dealing w/NAs
y4 <- c(svyby(~country_wrongtrack, ~ stay_leave, dsgn.stay19.tr, svymean)$country_wrongtrack,svyby(~country_wrongtrack, ~ stay_leave, dsgn.never19.tr, svymean)$country_wrongtrack)
y4_error <- c(sqrt(svyby(~country_wrongtrack, ~stay_leave, dsgn.stay19.tr, svyvar)$country_wrongtrack)/sqrt(nrow(dsgn.stay19.tr)),sqrt(svyby(~country_wrongtrack, ~stay_leave, dsgn.never19.tr, svyvar)$country_wrongtrack)/sqrt(nrow(dsgn.never19.tr)))
# trump_disapp_econ
y5 <- c(svyby(~trump_disapp_econ, ~ stay_leave, dsgn.stay19, svymean)$trump_disapp_econ,svyby(~trump_disapp_econ, ~ stay_leave, dsgn.never19, svymean)$trump_disapp_econ)
y5_error <- c(sqrt(svyby(~trump_disapp_econ, ~stay_leave, dsgn.stay19, svyvar)$trump_disapp_econ)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~trump_disapp_econ, ~stay_leave, dsgn.never19, svyvar)$trump_disapp_econ)/sqrt(nrow(d19_Jo_Nv)))
# trump_disapp_stimjobs 
y6 <- c(svyby(~trump_disapp_stimjobs, ~ stay_leave, dsgn.stay19, svymean)$trump_disapp_stimjobs,svyby(~trump_disapp_stimjobs, ~ stay_leave, dsgn.never19, svymean)$trump_disapp_stimjobs)
y6_error <- c(sqrt(svyby(~trump_disapp_stimjobs, ~stay_leave, dsgn.stay19, svyvar)$trump_disapp_stimjobs)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~trump_disapp_stimjobs, ~stay_leave, dsgn.never19, svyvar)$trump_disapp_stimjobs)/sqrt(nrow(d19_Jo_Nv)))
## combining calculations
df2 <- rbind.data.frame(cbind(x,y=y2,y_error=y2_error), cbind(x,y=y3,y_error=y3_error), cbind(x,y=y4,y_error=y4_error), cbind(x,y=y5,y_error=y5_error), cbind(x,y=y6,y_error=y6_error))
df2$y <- as.numeric(as.character(df2$y))
df2$y_error <- as.numeric(as.character(df2$y_error))
df2$measure <- c(rep("US economy is on wrong track",4),rep("US economy is very weak",4),rep("Country is on wrong track",4),rep("Disapprove of Trump on economy",4),rep("Disapprove of Trump on stimulating jobs",4))
colnames(df2) <- c("model","estimate","std.error","term")
# getting ordering by R Loyal - D Loyal difference
if(nrow(d19.sub)==907){
  diffs <- data.frame(matrix(ncol = 2, nrow = length(unique(df2$term))))
  names(diffs) <- c("term","diff")
  diffs$term <- unique(df2$term)
  for(i in 1:length(unique(df2$term))){
    subset <- subset(df2, term==unique(df2$term)[i])
    diffs[i,2] <- subset$estimate[subset$model=="Rep Loyalists"] - subset$estimate[subset$model=="Dem Loyalists"]
  }
  diffs$rank <- rank(abs(diffs$diff))
  diffs.sub.ecsoc <- subset(diffs, select = c("term","rank"))}
df2 <- merge(df2, diffs.sub.ecsoc, by = "term")
df2 <- rbind.data.frame(df2[df2$rank==5,],df2[df2$rank==4,],df2[df2$rank==3,],df2[df2$rank==2,],df2[df2$rank==1,])
## plotting data
dwplot(df2, dodge_size = 0.4, dot_args=list(aes(colour = model, shape = model), size = 3), whisker_args=list(aes(colour = model))) + 
  theme_minimal() + xlab("Sociotropic economic anxiety (least to most)") + 
  theme(legend.title = element_blank(), legend.position = "top",legend.key.width = unit(.5,"cm"), 
        plot.title = element_text(size = 12), axis.title.x = element_text(size = 10), axis.text = element_text(size=10)) + 
  xlim(-1,1) + scale_colour_manual(breaks=c("Dem Loyalists","Dem Switchers","Rep Switchers","Rep Loyalists"), values= c("darkorange", alpha("darkorange", .45), alpha("purple3", .45), "purple3")) +
  scale_shape_manual(breaks=c("Dem Loyalists","Dem Switchers","Rep Switchers","Rep Loyalists"), values = c(20,20,20,20))
### run only one of the next two lines!
ggsave("online_appendix_A/plots/FigF_1a_DEMS_econanxiety_socio.png", width = 7.8, height = 4.25)
#ggsave("online_appendix_A/plots/FigF_1b_REPS_econanxiety_socio.png", width = 7.8, height = 4.25)

#### FIGURE F - part 2 ####
## doing calculations for plots
# ECON_LADDER_past_future 
y2 <- c(svyby(~ECON_LADDER_past_future, ~ stay_leave, dsgn.stay19, svymean)$ECON_LADDER_past_future,svyby(~ECON_LADDER_past_future, ~ stay_leave, dsgn.never19, svymean)$ECON_LADDER_past_future)
y2_error <- c(sqrt(svyby(~ECON_LADDER_past_future, ~stay_leave, dsgn.stay19, svyvar)$ECON_LADDER_past_future)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~ECON_LADDER_past_future, ~stay_leave, dsgn.never19, svyvar)$ECON_LADDER_past_future)/sqrt(nrow(d19_Jo_Nv)))
# ECON_LADDER_past_present
y3 <- c(svyby(~ECON_LADDER_past_present, ~ stay_leave, dsgn.stay19, svymean)$ECON_LADDER_past_present,svyby(~ECON_LADDER_past_present, ~ stay_leave, dsgn.never19, svymean)$ECON_LADDER_past_present)
y3_error <- c(sqrt(svyby(~ECON_LADDER_past_present, ~stay_leave, dsgn.stay19, svyvar)$ECON_LADDER_past_present)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~ECON_LADDER_past_present, ~stay_leave, dsgn.never19, svyvar)$ECON_LADDER_past_present)/sqrt(nrow(d19_Jo_Nv)))
# ECON_LADDER_present_future
y4 <- c(svyby(~ECON_LADDER_present_future, ~ stay_leave, dsgn.stay19, svymean)$ECON_LADDER_present_future,svyby(~ECON_LADDER_present_future, ~ stay_leave, dsgn.never19, svymean)$ECON_LADDER_present_future)
y4_error <- c(sqrt(svyby(~ECON_LADDER_present_future, ~stay_leave, dsgn.stay19, svyvar)$ECON_LADDER_present_future)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~ECON_LADDER_present_future, ~stay_leave, dsgn.never19, svyvar)$ECON_LADDER_present_future)/sqrt(nrow(d19_Jo_Nv)))
# personalfinance
dsgn.stay19.pf <- subset(dsgn.stay19, !is.na(dsgn.stay19$variables$personalfinance)) # dealing w/NAs
dsgn.never19.pf <- subset(dsgn.never19, !is.na(dsgn.never19$variables$personalfinance)) # dealing w/NAs
y5 <- c(svyby(~personalfinance, ~ stay_leave, dsgn.stay19.pf, svymean)$personalfinance,svyby(~personalfinance, ~ stay_leave, dsgn.never19.pf, svymean)$personalfinance)
y5_error <- c(sqrt(svyby(~personalfinance, ~stay_leave, dsgn.stay19.pf, svyvar)$personalfinance)/sqrt(nrow(dsgn.stay19.pf)),sqrt(svyby(~personalfinance, ~stay_leave, dsgn.never19.pf, svyvar)$personalfinance)/sqrt(nrow(dsgn.never19.pf)))
y5 <- y5*-1 # reversing scale
# PARENTS_rec
dsgn.stay19.par <- subset(dsgn.stay19, !is.na(dsgn.stay19$variables$PARENTS_rec)) # dealing w/NAs
dsgn.never19.par <- subset(dsgn.never19, !is.na(dsgn.never19$variables$PARENTS_rec)) # dealing w/NAs
y6 <- c(svyby(~PARENTS_rec, ~ stay_leave, dsgn.stay19.par, svymean)$PARENTS_rec,svyby(~PARENTS_rec, ~ stay_leave, dsgn.never19.par, svymean)$PARENTS_rec)
y6_error <- c(sqrt(svyby(~PARENTS_rec, ~stay_leave, dsgn.stay19.par, svyvar)$PARENTS_rec)/sqrt(nrow(dsgn.stay19.par)),sqrt(svyby(~PARENTS_rec, ~stay_leave, dsgn.never19.par, svyvar)$PARENTS_rec)/sqrt(nrow(dsgn.never19.par)))
y6 <- y6*-1 # reversing scale
# CHILDREN_rec
dsgn.stay19.ch <- subset(dsgn.stay19, !is.na(dsgn.stay19$variables$CHILDREN_rec)) # dealing w/NAs
dsgn.never19.ch <- subset(dsgn.never19, !is.na(dsgn.never19$variables$CHILDREN_rec)) # dealing w/NAs
y7 <- c(svyby(~CHILDREN_rec, ~ stay_leave, dsgn.stay19.ch, svymean)$CHILDREN_rec,svyby(~CHILDREN_rec, ~ stay_leave, dsgn.never19.ch, svymean)$CHILDREN_rec)
y7_error <- c(sqrt(svyby(~CHILDREN_rec, ~stay_leave, dsgn.stay19.ch, svyvar)$CHILDREN_rec)/sqrt(nrow(dsgn.stay19.ch)),sqrt(svyby(~CHILDREN_rec, ~stay_leave, dsgn.never19.ch, svyvar)$CHILDREN_rec)/sqrt(nrow(dsgn.never19.ch)))
y7 <- y7*-1 # reversing scale
## combining calculations
df2 <- rbind.data.frame(cbind(x,y=y2,y_error=y2_error), cbind(x,y=y3,y_error=y3_error), cbind(x,y=y4,y_error=y4_error), cbind(x,y=y5,y_error=y5_error), cbind(x,y=y6,y_error=y6_error), cbind(x,y=y7,y_error=y7_error))
df2$y <- as.numeric(as.character(df2$y))
df2$y_error <- as.numeric(as.character(df2$y_error))
df2$measure <- c(rep("Family's finances as teen minus\nexpected of next gen",4),rep("Family's finances as teen\nminus current",4),rep("Current family finances minus\nexpected of next gen",4),rep("Personal financial situation is\ngetting worse",4),rep("Job opportunities for your children\nwill be worse",4),rep("Job opportunities for you worse\nthan parents'",4))
colnames(df2) <- c("model","estimate","std.error","term")
# getting ordering by R Loyal - D Loyal difference
if(nrow(d19.sub)==907){
  diffs <- data.frame(matrix(ncol = 2, nrow = length(unique(df2$term))))
  names(diffs) <- c("term","diff")
  diffs$term <- unique(df2$term)
  for(i in 1:length(unique(df2$term))){
    subset <- subset(df2, term==unique(df2$term)[i])
    diffs[i,2] <- subset$estimate[subset$model=="Rep Loyalists"] - subset$estimate[subset$model=="Dem Loyalists"]
  }
  diffs$rank <- rank(abs(diffs$diff))
  diffs.sub.ecself <- subset(diffs, select = c("term","rank"))}
df2 <- merge(df2, diffs.sub.ecself, by = "term")
df2 <- rbind.data.frame(df2[df2$rank==6,],df2[df2$rank==5,],df2[df2$rank==4,],df2[df2$rank==3,],df2[df2$rank==2,],df2[df2$rank==1,])
## plotting data
dwplot(df2, dodge_size = 0.4, dot_args=list(aes(colour = model, shape = model), size = 3), whisker_args=list(aes(colour = model))) + 
  theme_minimal() + xlab("Personal economic anxiety (least to most)") + 
  theme(legend.title = element_blank(), legend.position = "top",legend.key.width = unit(.5,"cm"), 
        plot.title = element_text(size = 12), axis.title.x = element_text(size = 10), axis.text = element_text(size=10)) + 
  xlim(-1,1) + scale_colour_manual(breaks=c("Dem Loyalists","Dem Switchers","Rep Switchers","Rep Loyalists"), values= c("darkorange", alpha("darkorange", .45), alpha("purple3", .45), "purple3")) +
  scale_shape_manual(breaks=c("Dem Loyalists","Dem Switchers","Rep Switchers","Rep Loyalists"), values = c(20,20,20,20))
### run only one of the next two lines!
ggsave("online_appendix_A/plots/FigF_2a_DEMS_econanxiety_self.png", width = 7.8, height = 4.6)
#ggsave("online_appendix_A/plots/FigF_2b_REPS_econanxiety_self.png", width = 7.8, height = 4.6)

#### FIGURE F - part 3 ####
## doing calculations for plots
# populism / taxes_scale_v2 
y5 <- c(svyby(~taxes_scale_v2, ~ stay_leave, dsgn.stay19, svymean)$taxes_scale_v2,svyby(~taxes_scale_v2, ~ stay_leave, dsgn.never19, svymean)$taxes_scale_v2)
y5_error <- c(sqrt(svyby(~taxes_scale_v2, ~stay_leave, dsgn.stay19, svyvar)$taxes_scale_v2)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~taxes_scale_v2, ~stay_leave, dsgn.never19, svyvar)$taxes_scale_v2)/sqrt(nrow(d19_Jo_Nv)))
y5 <- y5*-1 # reversing scale
# authoritarianism / BREAK_RULES_rec
y2 <- c(svyby(~BREAK_RULES_rec, ~ stay_leave, dsgn.stay19, svymean)$BREAK_RULES_rec,svyby(~BREAK_RULES_rec, ~ stay_leave, dsgn.never19, svymean)$BREAK_RULES_rec)
y2_error <- c(sqrt(svyby(~BREAK_RULES_rec, ~stay_leave, dsgn.stay19, svyvar)$BREAK_RULES_rec)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~BREAK_RULES_rec, ~stay_leave, dsgn.never19, svyvar)$BREAK_RULES_rec)/sqrt(nrow(d19_Jo_Nv)))
# COMPROMISE
y4 <- c(svyby(~COMPROMISE, ~ stay_leave, dsgn.stay19, svymean)$COMPROMISE,svyby(~COMPROMISE, ~ stay_leave, dsgn.never19, svymean)$COMPROMISE)
y4_error <- c(sqrt(svyby(~COMPROMISE, ~stay_leave, dsgn.stay19, svyvar)$COMPROMISE)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~COMPROMISE, ~stay_leave, dsgn.never19, svyvar)$COMPROMISE)/sqrt(nrow(d19_Jo_Nv)))
# Law favorability scale : wealthy people, large corporations
y6 <- c(svyby(~laws_toofavorable_wealthycorps, ~ stay_leave, dsgn.stay19, svymean)$laws_toofavorable_wealthycorps,svyby(~laws_toofavorable_wealthycorps, ~ stay_leave, dsgn.never19, svymean)$laws_toofavorable_wealthycorps)
y6 <- y6*-1 # reversing scale
y6_error <- c(sqrt(svyby(~laws_toofavorable_wealthycorps, ~stay_leave, dsgn.stay19, svyvar)$laws_toofavorable_wealthycorps)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~laws_toofavorable_wealthycorps, ~stay_leave, dsgn.never19, svyvar)$laws_toofavorable_wealthycorps)/sqrt(nrow(d19_Jo_Nv)))
## combining calculations
df2 <- rbind.data.frame(cbind(x,y=y6,y_error=y6_error),cbind(x,y=y5,y_error=y5_error), cbind(x,y=y2,y_error=y2_error), cbind(x,y=y4,y_error=y4_error))
df2$y <- as.numeric(as.character(df2$y))
df2$y_error <- as.numeric(as.character(df2$y_error))
df2$measure <- c(rep("Laws not favorable enough toward\nwealthy people, large corps",4),rep("Wealthy people, large corps receive less\nthan they contribute in taxes",4),rep("We need a leader willing to break rules",4),rep("More important for leader to stick to\nbeliefs than compromise",4))
colnames(df2) <- c("model","estimate","std.error","term")
# getting ordering by R Loyal - D Loyal difference
if(nrow(d19.sub)==907){
  diffs <- data.frame(matrix(ncol = 2, nrow = length(unique(df2$term))))
  names(diffs) <- c("term","diff")
  diffs$term <- unique(df2$term)
  for(i in 1:length(unique(df2$term))){
    subset <- subset(df2, term==unique(df2$term)[i])
    diffs[i,2] <- subset$estimate[subset$model=="Rep Loyalists"] - subset$estimate[subset$model=="Dem Loyalists"]
  }
  diffs$rank <- rank(diffs$diff)
  diffs.sub.gov <- subset(diffs, select = c("term","rank"))}
df2 <- merge(df2, diffs.sub.gov, by = "term")
df2 <- rbind.data.frame(df2[df2$rank==4,],df2[df2$rank==3,],df2[df2$rank==2,],df2[df2$rank==1,])
## plotting data
dwplot(df2, dodge_size = 0.4, dot_args=list(aes(colour = model, shape = model), size = 3), whisker_args=list(aes(colour = model))) + 
  theme_minimal() + xlab("Disagree to Agree") + 
  theme(legend.title = element_blank(), legend.position = "top",legend.key.width = unit(.5,"cm"), 
        plot.title = element_text(size = 12), axis.title.x = element_text(size = 10), axis.text = element_text(size=10)) + 
  xlim(-1,1) + scale_colour_manual(breaks=c("Dem Loyalists","Dem Switchers","Rep Switchers","Rep Loyalists"), values= c("darkorange", alpha("darkorange", .45), alpha("purple3", .45), "purple3")) +
  scale_shape_manual(breaks=c("Dem Loyalists","Dem Switchers","Rep Switchers","Rep Loyalists"), values = c(20,20,20,20))
### run only one of the next two lines!
ggsave("online_appendix_A/plots/FigF_3a_DEMS_governance_styles_new.png", width = 7.8, height = 3.3)
#ggsave("online_appendix_A/plots/FigF_3b_REPS_governance_styles_new.png", width = 7.8, height = 3.3)


#### FIGURE F - part 5 ####
## setting up datasets for is it racist? and racism/nativism plots
# subsetting datasets and standardizing scales
vars <- c("RESPID","weight","stay_leave","ptystr","RACISM_OKR1","RACISM_OKR2","RACISM_OKR3","RACISM_OKR4",
          "RACISM_OKR5","RACISM_OKR6","RACISM_OKR7","RACISM_OKR8","RACISM_OKR9","RACISM_OKR10",
          "is_it_racist_scale1","is_it_racist_scale2")
d17.sub <- d17[vars]
x <- data.frame(sapply(c(5:16), function(i){scales::rescale(as.numeric(unlist(d17.sub[,i])), to = c(-1,1))}))
d17.sub <- cbind.data.frame(d17.sub[,1:4],x)
colnames(d17.sub) <- vars
# subsetting among partisans
d17rep <- subset(d17.sub, ptystr<4)
d17dem <- subset(d17.sub, ptystr>4)

#### RUN THE FOLLOWING CODE FOR DEMOCRAT SUBSET FIRST, THEN REPUBLICAN 
#### TO REPLICATE FIGURE E, parts 4-5
### this way, the plotted items appear in the same order for Dem and Rep plots
### run only one of the next two lines!
d17.sub <- d17dem
#d17.sub <- d17rep

# subset datasets
d17_St_Lv <- subset(d17.sub, stay_leave==1 | stay_leave==2)
d17_Jo_Nv <- subset(d17.sub, stay_leave==3 | stay_leave==4)
# setting up survey designs
dsgn.stay <- svydesign(ids = ~RESPID, weights = ~weight, data = d17_St_Lv)
dsgn.never <- svydesign(ids = ~RESPID, weights = ~weight, data = d17_Jo_Nv)
# prep for plot DFs
x <- c("Rep Loyalists","Rep Switchers","Dem Switchers","Dem Loyalists")

## doing calculations for is it racist? plot
# RACISM_OKR1 
y1 <- c(svyby(~RACISM_OKR1, ~ stay_leave, dsgn.stay, svymean)$RACISM_OKR1,svyby(~RACISM_OKR1, ~ stay_leave, dsgn.never, svymean)$RACISM_OKR1)
y1_error <- c(sqrt(svyby(~RACISM_OKR1, ~stay_leave, dsgn.stay, svyvar)$RACISM_OKR1)/sqrt(nrow(d17_St_Lv)),sqrt(svyby(~RACISM_OKR1, ~stay_leave, dsgn.never, svyvar)$RACISM_OKR1)/sqrt(nrow(d17_Jo_Nv)))
# RACISM_OKR2
y2 <- c(svyby(~RACISM_OKR2, ~ stay_leave, dsgn.stay, svymean)$RACISM_OKR2,svyby(~RACISM_OKR2, ~ stay_leave, dsgn.never, svymean)$RACISM_OKR2)
y2_error <- c(sqrt(svyby(~RACISM_OKR2, ~stay_leave, dsgn.stay, svyvar)$RACISM_OKR2)/sqrt(nrow(d17_St_Lv)),sqrt(svyby(~RACISM_OKR2, ~stay_leave, dsgn.never, svyvar)$RACISM_OKR2)/sqrt(nrow(d17_Jo_Nv)))
# RACISM_OKR3
y3 <- c(svyby(~RACISM_OKR3, ~ stay_leave, dsgn.stay, svymean)$RACISM_OKR3,svyby(~RACISM_OKR3, ~ stay_leave, dsgn.never, svymean)$RACISM_OKR3)
y3_error <- c(sqrt(svyby(~RACISM_OKR3, ~stay_leave, dsgn.stay, svyvar)$RACISM_OKR3)/sqrt(nrow(d17_St_Lv)),sqrt(svyby(~RACISM_OKR3, ~stay_leave, dsgn.never, svyvar)$RACISM_OKR3)/sqrt(nrow(d17_Jo_Nv)))
# RACISM_OKR4
y4 <- c(svyby(~RACISM_OKR4, ~ stay_leave, dsgn.stay, svymean)$RACISM_OKR4,svyby(~RACISM_OKR4, ~ stay_leave, dsgn.never, svymean)$RACISM_OKR4)
y4_error <- c(sqrt(svyby(~RACISM_OKR4, ~stay_leave, dsgn.stay, svyvar)$RACISM_OKR4)/sqrt(nrow(d17_St_Lv)),sqrt(svyby(~RACISM_OKR4, ~stay_leave, dsgn.never, svyvar)$RACISM_OKR4)/sqrt(nrow(d17_Jo_Nv)))
# RACISM_OKR5
y5 <- c(svyby(~RACISM_OKR5, ~ stay_leave, dsgn.stay, svymean)$RACISM_OKR5,svyby(~RACISM_OKR5, ~ stay_leave, dsgn.never, svymean)$RACISM_OKR5)
y5_error <- c(sqrt(svyby(~RACISM_OKR5, ~stay_leave, dsgn.stay, svyvar)$RACISM_OKR5)/sqrt(nrow(d17_St_Lv)),sqrt(svyby(~RACISM_OKR5, ~stay_leave, dsgn.never, svyvar)$RACISM_OKR5)/sqrt(nrow(d17_Jo_Nv)))
# RACISM_OKR6
y6 <- c(svyby(~RACISM_OKR6, ~ stay_leave, dsgn.stay, svymean)$RACISM_OKR6,svyby(~RACISM_OKR6, ~ stay_leave, dsgn.never, svymean)$RACISM_OKR6)
y6_error <- c(sqrt(svyby(~RACISM_OKR6, ~stay_leave, dsgn.stay, svyvar)$RACISM_OKR6)/sqrt(nrow(d17_St_Lv)),sqrt(svyby(~RACISM_OKR6, ~stay_leave, dsgn.never, svyvar)$RACISM_OKR6)/sqrt(nrow(d17_Jo_Nv)))
# RACISM_OKR7 
y7 <- c(svyby(~RACISM_OKR7, ~ stay_leave, dsgn.stay, svymean)$RACISM_OKR7,svyby(~RACISM_OKR7, ~ stay_leave, dsgn.never, svymean)$RACISM_OKR7)
y7_error <- c(sqrt(svyby(~RACISM_OKR7, ~stay_leave, dsgn.stay, svyvar)$RACISM_OKR7)/sqrt(nrow(d17_St_Lv)),sqrt(svyby(~RACISM_OKR7, ~stay_leave, dsgn.never, svyvar)$RACISM_OKR7)/sqrt(nrow(d17_Jo_Nv)))
# RACISM_OKR8 
y8 <- c(svyby(~RACISM_OKR8, ~ stay_leave, dsgn.stay, svymean)$RACISM_OKR8,svyby(~RACISM_OKR8, ~ stay_leave, dsgn.never, svymean)$RACISM_OKR8)
y8_error <- c(sqrt(svyby(~RACISM_OKR8, ~stay_leave, dsgn.stay, svyvar)$RACISM_OKR8)/sqrt(nrow(d17_St_Lv)),sqrt(svyby(~RACISM_OKR8, ~stay_leave, dsgn.never, svyvar)$RACISM_OKR8)/sqrt(nrow(d17_Jo_Nv)))
# RACISM_OKR9 
y9 <- c(svyby(~RACISM_OKR9, ~ stay_leave, dsgn.stay, svymean)$RACISM_OKR9,svyby(~RACISM_OKR9, ~ stay_leave, dsgn.never, svymean)$RACISM_OKR9)
y9_error <- c(sqrt(svyby(~RACISM_OKR9, ~stay_leave, dsgn.stay, svyvar)$RACISM_OKR9)/sqrt(nrow(d17_St_Lv)),sqrt(svyby(~RACISM_OKR9, ~stay_leave, dsgn.never, svyvar)$RACISM_OKR9)/sqrt(nrow(d17_Jo_Nv)))
# RACISM_OKR10
y10 <- c(svyby(~RACISM_OKR10, ~ stay_leave, dsgn.stay, svymean)$RACISM_OKR10,svyby(~RACISM_OKR10, ~ stay_leave, dsgn.never, svymean)$RACISM_OKR10)
y10_error <- c(sqrt(svyby(~RACISM_OKR10, ~stay_leave, dsgn.stay, svyvar)$RACISM_OKR10)/sqrt(nrow(d17_St_Lv)),sqrt(svyby(~RACISM_OKR10, ~stay_leave, dsgn.never, svyvar)$RACISM_OKR10)/sqrt(nrow(d17_Jo_Nv)))
## combining calculations 
df2 <- rbind.data.frame(cbind(x,y=y4,y_error=y4_error),cbind(x,y=y2,y_error=y2_error),cbind(x,y=y7,y_error=y7_error),cbind(x,y=y3,y_error=y3_error),cbind(x,y=y9,y_error=y9_error),cbind(x,y=y6,y_error=y6_error),cbind(x,y=y1,y_error=y1_error),cbind(x,y=y5,y_error=y5_error),cbind(x,y=y8,y_error=y8_error),cbind(x,y=y10,y_error=y10_error))
df2$y <- as.numeric(as.character(df2$y))
df2$y_error <- as.numeric(as.character(df2$y_error))
df2$measure <- c(rep("Voting for Donald Trump",4),rep("Saying immigrants commit too many crimes",4),
                 rep("Wanting to fly the Confederate flag",4),rep("Saying people in America should\nspeak English",4),
                 rep("Agreeing when police shoot civilians,\nit's for a good reason",4),rep("Using a word some see as offensive\nabout racial group",4),
                 rep("Believing country was a better place to live\nin the past",4),rep("Telling a joke about a racial group",4),
                 rep("Wanting your children to attend school\nwith others from own cultural background",4),rep("Agreeing welfare recipients should work\nto get benefits",4))
colnames(df2) <- c("model","estimate","std.error","term")
# getting ordering by R Loyal - D Loyal difference
if(nrow(d17.sub)==1004){
  diffs <- data.frame(matrix(ncol = 2, nrow = length(unique(df2$term))))
  names(diffs) <- c("term","diff")
  diffs$term <- unique(df2$term)
  for(i in 1:length(unique(df2$term))){
    subset <- subset(df2, term==unique(df2$term)[i])
    diffs[i,2] <- subset$estimate[subset$model=="Rep Loyalists"] - subset$estimate[subset$model=="Dem Loyalists"]
  }
  diffs$rank <- rank(abs(diffs$diff))
  diffs.sub.race <- subset(diffs, select = c("term","rank"))}
df2 <- merge(df2, diffs.sub.race, by = "term")
df2 <- rbind.data.frame(df2[df2$rank==10,],df2[df2$rank==9,],df2[df2$rank==8,],df2[df2$rank==7,],df2[df2$rank==6,],df2[df2$rank==5,],df2[df2$rank==4,],df2[df2$rank==3,],df2[df2$rank==2,],df2[df2$rank==1,])
## plotting data
dwplot(df2, dodge_size = 0.4, dot_args=list(aes(colour = model, shape = model), size = 3), whisker_args=list(aes(colour = model))) + 
  theme_minimal() + xlab("Agree statement is racist -to- Agree statement is not racist") + 
  theme(legend.title = element_blank(), legend.position = "top",legend.key.width = unit(.5,"cm"), 
        plot.title = element_text(size = 12), axis.title.x = element_text(size = 10), axis.text = element_text(size=10)) + 
  xlim(-1,1) + scale_colour_manual(breaks=c("Dem Loyalists","Dem Switchers","Rep Switchers","Rep Loyalists"), values= c("darkorange", alpha("darkorange", .45), alpha("purple3", .45), "purple3")) +
  scale_shape_manual(breaks=c("Dem Loyalists","Dem Switchers","Rep Switchers","Rep Loyalists"), values = c(20,20,20,20))
### run only one of the next two lines!
ggsave("online_appendix_A/plots/FigF_5a_DEMS_is_it_racist.png", width = 7.8, height = 5.9)
#ggsave("online_appendix_A/plots/FigF_5b_REPS_is_it_racist.png", width = 7.8, height = 5.9)


#### FIGURE F - part F ####
## doing calculations for plots
# cosmopolitanism scale
y2 <- c(svyby(~cosmo_scale2, ~ stay_leave, dsgn.stay19, svymean)$cosmo_scale2,svyby(~cosmo_scale2, ~ stay_leave, dsgn.never19, svymean)$cosmo_scale2)
y2 <- y2*-1 # reversing scale
y2_error <- c(sqrt(svyby(~cosmo_scale2, ~stay_leave, dsgn.stay19, svyvar)$cosmo_scale2)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~cosmo_scale2, ~stay_leave, dsgn.never19, svyvar)$cosmo_scale2)/sqrt(nrow(d19_Jo_Nv)))
# Truly American scale v1
y3 <- c(svyby(~american_scale_v1, ~ stay_leave, dsgn.stay19, svymean)$american_scale_v1,svyby(~american_scale_v1, ~ stay_leave, dsgn.never19, svymean)$american_scale_v1)
y3_error <- c(sqrt(svyby(~american_scale_v1, ~stay_leave, dsgn.stay19, svyvar)$american_scale_v1)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~american_scale_v1, ~stay_leave, dsgn.never19, svyvar)$american_scale_v1)/sqrt(nrow(d19_Jo_Nv)))
# Truly American scale v2 
y4 <- c(svyby(~american_scale_v2, ~ stay_leave, dsgn.stay19, svymean)$american_scale_v2,svyby(~american_scale_v2, ~ stay_leave, dsgn.never19, svymean)$american_scale_v2)
y4_error <- c(sqrt(svyby(~american_scale_v2, ~stay_leave, dsgn.stay19, svyvar)$american_scale_v2)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~american_scale_v2, ~stay_leave, dsgn.never19, svyvar)$american_scale_v2)/sqrt(nrow(d19_Jo_Nv)))
# Taxes scale v1 
y5 <- c(svyby(~taxes_scale_v1, ~ stay_leave, dsgn.stay19, svymean)$taxes_scale_v1,svyby(~taxes_scale_v1, ~ stay_leave, dsgn.never19, svymean)$taxes_scale_v1)
y5_error <- c(sqrt(svyby(~taxes_scale_v1, ~stay_leave, dsgn.stay19, svyvar)$taxes_scale_v1)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~taxes_scale_v1, ~stay_leave, dsgn.never19, svyvar)$taxes_scale_v1)/sqrt(nrow(d19_Jo_Nv)))
# Law favorability scale : blacks, hispanics, illegal immigrants 
y6 <- c(svyby(~laws_toofavorable, ~ stay_leave, dsgn.stay19, svymean)$laws_toofavorable,svyby(~laws_toofavorable, ~ stay_leave, dsgn.never19, svymean)$laws_toofavorable)
y6_error <- c(sqrt(svyby(~laws_toofavorable, ~stay_leave, dsgn.stay19, svyvar)$laws_toofavorable)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~laws_toofavorable, ~stay_leave, dsgn.never19, svyvar)$laws_toofavorable)/sqrt(nrow(d19_Jo_Nv)))
# Approval on Trump's job w/immigration / trump_immig_app
y7 <- c(svyby(~trump_immig_app, ~ stay_leave, dsgn.stay19, svymean)$trump_immig_app,svyby(~trump_immig_app, ~ stay_leave, dsgn.never19, svymean)$trump_immig_app)
y7 <- y7*-1 # reversing scale 
y7_error <- c(sqrt(svyby(~trump_immig_app, ~stay_leave, dsgn.stay19, svyvar)$trump_immig_app)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~trump_immig_app, ~stay_leave, dsgn.never19, svyvar)$trump_immig_app)/sqrt(nrow(d19_Jo_Nv)))
# racial resentment scale 
y8 <- c(svyby(~racial_resent_scale, ~ stay_leave, dsgn.stay19, svymean)$racial_resent_scale,svyby(~racial_resent_scale, ~ stay_leave, dsgn.never19, svymean)$racial_resent_scale)
y8_error <- c(sqrt(svyby(~racial_resent_scale, ~stay_leave, dsgn.stay19, svyvar)$racial_resent_scale)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~racial_resent_scale, ~stay_leave, dsgn.never19, svyvar)$racial_resent_scale)/sqrt(nrow(d19_Jo_Nv)))
# is it racist 1 - country better in past, people should speak English, voting for Trump, it's for good reason when police shoot civilians
y9 <- c(svyby(~is_it_racist_scale1, ~ stay_leave, dsgn.stay, svymean)$is_it_racist_scale1,svyby(~is_it_racist_scale1, ~ stay_leave, dsgn.never, svymean)$is_it_racist_scale1)
y9_error <- c(sqrt(svyby(~is_it_racist_scale1, ~stay_leave, dsgn.stay, svyvar)$is_it_racist_scale1)/sqrt(nrow(d17_St_Lv)),sqrt(svyby(~is_it_racist_scale1, ~stay_leave, dsgn.never, svyvar)$is_it_racist_scale1)/sqrt(nrow(d17_Jo_Nv)))
# is it racist 2 - telling joke about racial group, using offensive word for racial group
y10 <- c(svyby(~is_it_racist_scale2, ~ stay_leave, dsgn.stay, svymean)$is_it_racist_scale2,svyby(~is_it_racist_scale2, ~ stay_leave, dsgn.never, svymean)$is_it_racist_scale2)
y10_error <- c(sqrt(svyby(~is_it_racist_scale2, ~stay_leave, dsgn.stay, svyvar)$is_it_racist_scale2)/sqrt(nrow(d17_St_Lv)),sqrt(svyby(~is_it_racist_scale2, ~stay_leave, dsgn.never, svyvar)$is_it_racist_scale2)/sqrt(nrow(d17_Jo_Nv)))
## combining calculations
df2 <- rbind.data.frame(cbind(x,y=y2,y_error=y2_error), cbind(x,y=y3,y_error=y3_error), cbind(x,y=y4,y_error=y4_error), cbind(x,y=y5,y_error=y5_error), cbind(x,y=y6,y_error=y6_error), cbind(x,y=y7,y_error=y7_error),
                        cbind(x,y=y8,y_error=y8_error), cbind(x,y=y9,y_error=y9_error), cbind(x,y=y10,y_error=y10_error))
df2$y <- as.numeric(as.character(df2$y))
df2$y_error <- as.numeric(as.character(df2$y_error))
df2$measure <- c(rep("Cosmopolitanism",4),rep("Truly American if feel American,\nhave citizenship, respect institutions",4),rep("Truly American if born in US, lived in US\nmost of life, are Christian",4),
                 rep("Blacks, Hispanics, illegal immigrants receive\nmore than contribute in taxes",4),rep("Laws too favorable toward\nblacks, Hispanics, illegal immigrants",4),
                 rep("Disapproval of Trump handling immigration",4),rep("Racial resentment scale",4),rep("Not racist to say country better in past, people in\nU.S. should speak English; vote for Trump; agree\nwhen police shoot civilians, it's for good reason", 4),
                 rep("Not racist to tell joke about racial group,\nuse offensive word for racial group",4))
colnames(df2) <- c("model","estimate","std.error","term")
# getting ordering by R Loyal - D Loyal difference
if(nrow(d19.sub)==907){
  diffs <- data.frame(matrix(ncol = 2, nrow = length(unique(df2$term))))
  names(diffs) <- c("term","diff")
  diffs$term <- unique(df2$term)
  for(i in 1:length(unique(df2$term))){
    subset <- subset(df2, term==unique(df2$term)[i])
    diffs[i,2] <- subset$estimate[subset$model=="Rep Loyalists"] - subset$estimate[subset$model=="Dem Loyalists"]
  }
  diffs$rank <- rank(diffs$diff)
  diffs.sub.nat <- subset(diffs, select = c("term","rank"))}
df2 <- merge(df2, diffs.sub.nat, by = "term")
df2 <- rbind.data.frame(df2[df2$rank==9,],df2[df2$rank==8,],df2[df2$rank==7,],df2[df2$rank==6,],df2[df2$rank==5,],df2[df2$rank==4,],df2[df2$rank==3,],df2[df2$rank==2,],df2[df2$rank==1,])
## plotting data
dwplot(df2, dodge_size = 0.4, dot_args=list(aes(colour = model, shape = model), size = 3), whisker_args=list(aes(colour = model))) + 
  theme_minimal() + xlab("Left to Right") + 
  theme(legend.title = element_blank(), legend.position = "top",legend.key.width = unit(.5,"cm"), 
        plot.title = element_text(size = 12), axis.title.x = element_text(size = 10), axis.text = element_text(size=10)) + 
  xlim(-1,1) + scale_colour_manual(breaks=c("Dem Loyalists","Dem Switchers","Rep Switchers","Rep Loyalists"), values= c("darkorange", alpha("darkorange", .45), alpha("purple3", .45), "purple3")) +
  scale_shape_manual(breaks=c("Dem Loyalists","Dem Switchers","Rep Switchers","Rep Loyalists"), values = c(20,20,20,20))
### run only one of the next two lines!
ggsave("online_appendix_A/plots/FigF_4a_DEMS_nativism_and_racism.png", width = 8.5, height = 5)
#ggsave("online_appendix_A/plots/FigF_4b_REPS_nativism_and_racism.png", width = 8.5, height = 5)

#### GO BACK TO < RUN THE FOLLOWING CODE FOR DEMOCRAT SUBSET FIRST, THEN REPUBLICAN >


#### FIGURE G ####
## setting up datasets for plots
# subsetting datasets and standardizing scales
vars <- c("cosmo_scale2","american_scale_v1","american_scale_v2","taxes_scale_v1","laws_toofavorable","laws_toofavorable_wealthycorps", "trump_immig_app", #nativism
          "ECON_LADDER_past_future","ECON_LADDER_past_present","ECON_LADDER_present_future","personalfinance","PARENTS_rec","CHILDREN_rec", #economic anxiety (egocentric)
          "econ_wrongtrack","econ_strength", "country_wrongtrack","trump_disapp_econ","trump_disapp_stimjobs", #economic anxiety (sociotropic)
          "BREAK_RULES_rec","US_DARK_rec","COMPROMISE","taxes_scale_v2", #authoritarisnism
          "racial_resent_scale", # racism
          "TAXESR1","TAXESR2","TAXESR3","TAXESR4","TAXESR5","TAXESR6","TAXESR7","TAXESR8","TAXESR9","TAXESR10",
          "TAXESR11","TAXESR12","TAXESR13","TAXESR14","TAXESR15","LEGISLATIONR1","LEGISLATIONR2","LEGISLATIONR3",
          "LEGISLATIONR4","LEGISLATIONR5","LEGISLATIONR6","LEGISLATIONR7","LEGISLATIONR8","LEGISLATIONR9",             
          "LEGISLATIONR10","LEGISLATIONR11","LEGISLATIONR12","LEGISLATIONR13","LEGISLATIONR14","LEGISLATIONR15",
          "stay_leave","weight","RECORD",
          "ptystr")
d19.sub <- d19[vars]
x <- data.frame(sapply(c(1:53), function(i){scales::rescale(as.numeric(unlist(d19.sub[,i])), to = c(-1,1))}))
d19.sub <- cbind.data.frame(x,d19.sub[,54:57])
colnames(d19.sub) <- vars
# subset datasets
d19_St_Lv <- subset(d19.sub, stay_leave==1 | stay_leave==2)
d19_Jo_Nv <- subset(d19.sub, stay_leave==3 | stay_leave==4)
# setting up survey designs
dsgn.stay19 <- svydesign(ids = ~RECORD, weights = ~weight, data = d19_St_Lv)
dsgn.never19 <- svydesign(ids = ~RECORD, weights = ~weight, data = d19_Jo_Nv)
# prep for plot DFs
x <- c("Rep Loyalists","Rep Switchers","Dem Switchers","Dem Loyalists")

## doing calculations for plots, beginning w/taxes items
# TAXESR1
y <- c(svyby(~TAXESR1, ~ stay_leave, dsgn.stay19, svymean)$TAXESR1,svyby(~TAXESR1, ~ stay_leave, dsgn.never19, svymean)$TAXESR1)
y_error <- c(sqrt(svyby(~TAXESR1, ~stay_leave, dsgn.stay19, svyvar)$TAXESR1)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~TAXESR1, ~stay_leave, dsgn.never19, svyvar)$TAXESR1)/sqrt(nrow(d19_Jo_Nv)))
# TAXESR2 
y2 <- c(svyby(~TAXESR2, ~ stay_leave, dsgn.stay19, svymean)$TAXESR2,svyby(~TAXESR2, ~ stay_leave, dsgn.never19, svymean)$TAXESR2)
y2_error <- c(sqrt(svyby(~TAXESR2, ~stay_leave, dsgn.stay19, svyvar)$TAXESR2)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~TAXESR2, ~stay_leave, dsgn.never19, svyvar)$TAXESR2)/sqrt(nrow(d19_Jo_Nv)))
# TAXESR3
y3 <- c(svyby(~TAXESR3, ~ stay_leave, dsgn.stay19, svymean)$TAXESR3,svyby(~TAXESR3, ~ stay_leave, dsgn.never19, svymean)$TAXESR3)
y3_error <- c(sqrt(svyby(~TAXESR3, ~stay_leave, dsgn.stay19, svyvar)$TAXESR3)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~TAXESR3, ~stay_leave, dsgn.never19, svyvar)$TAXESR3)/sqrt(nrow(d19_Jo_Nv)))
# TAXESR4
y4 <- c(svyby(~TAXESR4, ~ stay_leave, dsgn.stay19, svymean)$TAXESR4,svyby(~TAXESR4, ~ stay_leave, dsgn.never19, svymean)$TAXESR4)
y4_error <- c(sqrt(svyby(~TAXESR4, ~stay_leave, dsgn.stay19, svyvar)$TAXESR4)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~TAXESR4, ~stay_leave, dsgn.never19, svyvar)$TAXESR4)/sqrt(nrow(d19_Jo_Nv)))
# TAXESR5
y5 <- c(svyby(~TAXESR5, ~ stay_leave, dsgn.stay19, svymean)$TAXESR5,svyby(~TAXESR5, ~ stay_leave, dsgn.never19, svymean)$TAXESR5)
y5_error <- c(sqrt(svyby(~TAXESR5, ~stay_leave, dsgn.stay19, svyvar)$TAXESR5)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~TAXESR5, ~stay_leave, dsgn.never19, svyvar)$TAXESR5)/sqrt(nrow(d19_Jo_Nv)))
# TAXESR6
y6 <- c(svyby(~TAXESR6, ~ stay_leave, dsgn.stay19, svymean)$TAXESR6,svyby(~TAXESR6, ~ stay_leave, dsgn.never19, svymean)$TAXESR6)
y6_error <- c(sqrt(svyby(~TAXESR6, ~stay_leave, dsgn.stay19, svyvar)$TAXESR6)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~TAXESR6, ~stay_leave, dsgn.never19, svyvar)$TAXESR6)/sqrt(nrow(d19_Jo_Nv)))
# TAXESR7
y7 <- c(svyby(~TAXESR7, ~ stay_leave, dsgn.stay19, svymean)$TAXESR7,svyby(~TAXESR7, ~ stay_leave, dsgn.never19, svymean)$TAXESR7)
y7_error <- c(sqrt(svyby(~TAXESR7, ~stay_leave, dsgn.stay19, svyvar)$TAXESR7)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~TAXESR7, ~stay_leave, dsgn.never19, svyvar)$TAXESR7)/sqrt(nrow(d19_Jo_Nv)))
# TAXESR8
y8 <- c(svyby(~TAXESR8, ~ stay_leave, dsgn.stay19, svymean)$TAXESR8,svyby(~TAXESR8, ~ stay_leave, dsgn.never19, svymean)$TAXESR8)
y8_error <- c(sqrt(svyby(~TAXESR8, ~stay_leave, dsgn.stay19, svyvar)$TAXESR8)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~TAXESR8, ~stay_leave, dsgn.never19, svyvar)$TAXESR8)/sqrt(nrow(d19_Jo_Nv)))
# TAXESR9
y9 <- c(svyby(~TAXESR9, ~ stay_leave, dsgn.stay19, svymean)$TAXESR9,svyby(~TAXESR9, ~ stay_leave, dsgn.never19, svymean)$TAXESR9)
y9_error <- c(sqrt(svyby(~TAXESR9, ~stay_leave, dsgn.stay19, svyvar)$TAXESR9)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~TAXESR9, ~stay_leave, dsgn.never19, svyvar)$TAXESR9)/sqrt(nrow(d19_Jo_Nv)))
# TAXESR10
y10 <- c(svyby(~TAXESR10, ~ stay_leave, dsgn.stay19, svymean)$TAXESR10,svyby(~TAXESR10, ~ stay_leave, dsgn.never19, svymean)$TAXESR10)
y10_error <- c(sqrt(svyby(~TAXESR10, ~stay_leave, dsgn.stay19, svyvar)$TAXESR10)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~TAXESR10, ~stay_leave, dsgn.never19, svyvar)$TAXESR10)/sqrt(nrow(d19_Jo_Nv)))
# TAXESR11
y11 <- c(svyby(~TAXESR11, ~ stay_leave, dsgn.stay19, svymean)$TAXESR11,svyby(~TAXESR11, ~ stay_leave, dsgn.never19, svymean)$TAXESR11)
y11_error <- c(sqrt(svyby(~TAXESR11, ~stay_leave, dsgn.stay19, svyvar)$TAXESR11)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~TAXESR11, ~stay_leave, dsgn.never19, svyvar)$TAXESR11)/sqrt(nrow(d19_Jo_Nv)))
# TAXESR12
y12 <- c(svyby(~TAXESR12, ~ stay_leave, dsgn.stay19, svymean)$TAXESR12,svyby(~TAXESR12, ~ stay_leave, dsgn.never19, svymean)$TAXESR12)
y12_error <- c(sqrt(svyby(~TAXESR12, ~stay_leave, dsgn.stay19, svyvar)$TAXESR12)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~TAXESR12, ~stay_leave, dsgn.never19, svyvar)$TAXESR12)/sqrt(nrow(d19_Jo_Nv)))
# TAXESR13
y13 <- c(svyby(~TAXESR13, ~ stay_leave, dsgn.stay19, svymean)$TAXESR13,svyby(~TAXESR13, ~ stay_leave, dsgn.never19, svymean)$TAXESR13)
y13_error <- c(sqrt(svyby(~TAXESR13, ~stay_leave, dsgn.stay19, svyvar)$TAXESR13)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~TAXESR13, ~stay_leave, dsgn.never19, svyvar)$TAXESR13)/sqrt(nrow(d19_Jo_Nv)))
# TAXESR14
y14 <- c(svyby(~TAXESR14, ~ stay_leave, dsgn.stay19, svymean)$TAXESR14,svyby(~TAXESR14, ~ stay_leave, dsgn.never19, svymean)$TAXESR14)
y14_error <- c(sqrt(svyby(~TAXESR14, ~stay_leave, dsgn.stay19, svyvar)$TAXESR14)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~TAXESR14, ~stay_leave, dsgn.never19, svyvar)$TAXESR14)/sqrt(nrow(d19_Jo_Nv)))
# TAXESR15
y15 <- c(svyby(~TAXESR15, ~ stay_leave, dsgn.stay19, svymean)$TAXESR15,svyby(~TAXESR15, ~ stay_leave, dsgn.never19, svymean)$TAXESR15)
y15_error <- c(sqrt(svyby(~TAXESR15, ~stay_leave, dsgn.stay19, svyvar)$TAXESR15)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~TAXESR15, ~stay_leave, dsgn.never19, svyvar)$TAXESR15)/sqrt(nrow(d19_Jo_Nv)))
## combining calculations
df2 <- rbind.data.frame(cbind(x,y=y,y_error=y_error),cbind(x,y=y2,y_error=y2_error),cbind(x,y=y3,y_error=y3_error),
                        cbind(x,y=y4,y_error=y4_error),cbind(x,y=y5,y_error=y5_error),cbind(x,y=y6,y_error=y6_error),
                        cbind(x,y=y7,y_error=y7_error),cbind(x,y=y8,y_error=y8_error),cbind(x,y=y9,y_error=y9_error),
                        cbind(x,y=y10,y_error=y10_error),cbind(x,y=y11,y_error=y11_error),cbind(x,y=y12,y_error=y12_error),
                        cbind(x,y=y13,y_error=y13_error),cbind(x,y=y14,y_error=y14_error),cbind(x,y=y15,y_error=y15_error))
df2$y <- -1*as.numeric(as.character(df2$y))
df2$y_error <- -1*as.numeric(as.character(df2$y_error))
df2$measure <- c(rep("Whites",4),rep("Blacks",4),rep("Hispanics",4),rep("People who live in large cities",4),
                 rep("People in small towns/rural",4),rep("Poor people",4),
                 rep("Middle class people",4),rep("Wealthy people",4),rep("Legal immigrants",4),
                 rep("Illegal immigrants",4),rep("Young adults",4),rep("The elderly",4),
                 rep("Large corporations",4),rep("Colleges and universities",4),rep("Churches/places of worship",4))
colnames(df2) <- c("model","estimate","std.error","term")
# getting ordering by R Loyal - D Loyal difference
diffs <- data.frame(matrix(ncol = 2, nrow = length(unique(df2$term))))
names(diffs) <- c("term","diff")
diffs$term <- unique(df2$term)
for(i in 1:length(unique(df2$term))){
  subset <- subset(df2, term==unique(df2$term)[i])
  diffs[i,2] <- subset$estimate[subset$model=="Rep Loyalists"] - subset$estimate[subset$model=="Dem Loyalists"]
}
diffs$rank <- rank(diffs$diff)
diffs.sub <- subset(diffs, select = c("term","rank"))
df2 <- merge(df2, diffs.sub, by = "term")
df2 <- rbind.data.frame(df2[df2$rank==15,],df2[df2$rank==14,],df2[df2$rank==13,],df2[df2$rank==12,],
                        df2[df2$rank==11,],df2[df2$rank==10,],df2[df2$rank==9,],df2[df2$rank==8,],
                        df2[df2$rank==7,],df2[df2$rank==6,],df2[df2$rank==5,],df2[df2$rank==4,],
                        df2[df2$rank==3,],df2[df2$rank==2,],df2[df2$rank==1,])
## doing calculations for plots, continuing w/legislation items
# LEGISLATIONR1
y <- c(svyby(~LEGISLATIONR1, ~ stay_leave, dsgn.stay19, svymean)$LEGISLATIONR1,svyby(~LEGISLATIONR1, ~ stay_leave, dsgn.never19, svymean)$LEGISLATIONR1)
y_error <- c(sqrt(svyby(~LEGISLATIONR1, ~stay_leave, dsgn.stay19, svyvar)$LEGISLATIONR1)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~LEGISLATIONR1, ~stay_leave, dsgn.never19, svyvar)$LEGISLATIONR1)/sqrt(nrow(d19_Jo_Nv)))
# LEGISLATIONR2
y2 <- c(svyby(~LEGISLATIONR2, ~ stay_leave, dsgn.stay19, svymean)$LEGISLATIONR2,svyby(~LEGISLATIONR2, ~ stay_leave, dsgn.never19, svymean)$LEGISLATIONR2)
y2_error <- c(sqrt(svyby(~LEGISLATIONR2, ~stay_leave, dsgn.stay19, svyvar)$LEGISLATIONR2)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~LEGISLATIONR2, ~stay_leave, dsgn.never19, svyvar)$LEGISLATIONR2)/sqrt(nrow(d19_Jo_Nv)))
# LEGISLATIONR3
y3 <- c(svyby(~LEGISLATIONR3, ~ stay_leave, dsgn.stay19, svymean)$LEGISLATIONR3,svyby(~LEGISLATIONR3, ~ stay_leave, dsgn.never19, svymean)$LEGISLATIONR3)
y3_error <- c(sqrt(svyby(~LEGISLATIONR3, ~stay_leave, dsgn.stay19, svyvar)$LEGISLATIONR3)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~LEGISLATIONR3, ~stay_leave, dsgn.never19, svyvar)$LEGISLATIONR3)/sqrt(nrow(d19_Jo_Nv)))
# LEGISLATIONR4
y4 <- c(svyby(~LEGISLATIONR4, ~ stay_leave, dsgn.stay19, svymean)$LEGISLATIONR4,svyby(~LEGISLATIONR4, ~ stay_leave, dsgn.never19, svymean)$LEGISLATIONR4)
y4_error <- c(sqrt(svyby(~LEGISLATIONR4, ~stay_leave, dsgn.stay19, svyvar)$LEGISLATIONR4)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~LEGISLATIONR4, ~stay_leave, dsgn.never19, svyvar)$LEGISLATIONR4)/sqrt(nrow(d19_Jo_Nv)))
# LEGISLATIONR5
y5 <- c(svyby(~LEGISLATIONR5, ~ stay_leave, dsgn.stay19, svymean)$LEGISLATIONR5,svyby(~LEGISLATIONR5, ~ stay_leave, dsgn.never19, svymean)$LEGISLATIONR5)
y5_error <- c(sqrt(svyby(~LEGISLATIONR5, ~stay_leave, dsgn.stay19, svyvar)$LEGISLATIONR5)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~LEGISLATIONR5, ~stay_leave, dsgn.never19, svyvar)$LEGISLATIONR5)/sqrt(nrow(d19_Jo_Nv)))
# LEGISLATIONR6
y6 <- c(svyby(~LEGISLATIONR6, ~ stay_leave, dsgn.stay19, svymean)$LEGISLATIONR6,svyby(~LEGISLATIONR6, ~ stay_leave, dsgn.never19, svymean)$LEGISLATIONR6)
y6_error <- c(sqrt(svyby(~LEGISLATIONR6, ~stay_leave, dsgn.stay19, svyvar)$LEGISLATIONR6)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~LEGISLATIONR6, ~stay_leave, dsgn.never19, svyvar)$LEGISLATIONR6)/sqrt(nrow(d19_Jo_Nv)))
# LEGISLATIONR7
y7 <- c(svyby(~LEGISLATIONR7, ~ stay_leave, dsgn.stay19, svymean)$LEGISLATIONR7,svyby(~LEGISLATIONR7, ~ stay_leave, dsgn.never19, svymean)$LEGISLATIONR7)
y7_error <- c(sqrt(svyby(~LEGISLATIONR7, ~stay_leave, dsgn.stay19, svyvar)$LEGISLATIONR7)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~LEGISLATIONR7, ~stay_leave, dsgn.never19, svyvar)$LEGISLATIONR7)/sqrt(nrow(d19_Jo_Nv)))
# LEGISLATIONR8
y8 <- c(svyby(~LEGISLATIONR8, ~ stay_leave, dsgn.stay19, svymean)$LEGISLATIONR8,svyby(~LEGISLATIONR8, ~ stay_leave, dsgn.never19, svymean)$LEGISLATIONR8)
y8_error <- c(sqrt(svyby(~LEGISLATIONR8, ~stay_leave, dsgn.stay19, svyvar)$LEGISLATIONR8)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~LEGISLATIONR8, ~stay_leave, dsgn.never19, svyvar)$LEGISLATIONR8)/sqrt(nrow(d19_Jo_Nv)))
# LEGISLATIONR9
y9 <- c(svyby(~LEGISLATIONR9, ~ stay_leave, dsgn.stay19, svymean)$LEGISLATIONR9,svyby(~LEGISLATIONR9, ~ stay_leave, dsgn.never19, svymean)$LEGISLATIONR9)
y9_error <- c(sqrt(svyby(~LEGISLATIONR9, ~stay_leave, dsgn.stay19, svyvar)$LEGISLATIONR9)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~LEGISLATIONR9, ~stay_leave, dsgn.never19, svyvar)$LEGISLATIONR9)/sqrt(nrow(d19_Jo_Nv)))
# LEGISLATIONR10
y10 <- c(svyby(~LEGISLATIONR10, ~ stay_leave, dsgn.stay19, svymean)$LEGISLATIONR10,svyby(~LEGISLATIONR10, ~ stay_leave, dsgn.never19, svymean)$LEGISLATIONR10)
y10_error <- c(sqrt(svyby(~LEGISLATIONR10, ~stay_leave, dsgn.stay19, svyvar)$LEGISLATIONR10)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~LEGISLATIONR10, ~stay_leave, dsgn.never19, svyvar)$LEGISLATIONR10)/sqrt(nrow(d19_Jo_Nv)))
# LEGISLATIONR11
y11 <- c(svyby(~LEGISLATIONR11, ~ stay_leave, dsgn.stay19, svymean)$LEGISLATIONR11,svyby(~LEGISLATIONR11, ~ stay_leave, dsgn.never19, svymean)$LEGISLATIONR11)
y11_error <- c(sqrt(svyby(~LEGISLATIONR11, ~stay_leave, dsgn.stay19, svyvar)$LEGISLATIONR11)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~LEGISLATIONR11, ~stay_leave, dsgn.never19, svyvar)$LEGISLATIONR11)/sqrt(nrow(d19_Jo_Nv)))
# LEGISLATIONR12
y12 <- c(svyby(~LEGISLATIONR12, ~ stay_leave, dsgn.stay19, svymean)$LEGISLATIONR12,svyby(~LEGISLATIONR12, ~ stay_leave, dsgn.never19, svymean)$LEGISLATIONR12)
y12_error <- c(sqrt(svyby(~LEGISLATIONR12, ~stay_leave, dsgn.stay19, svyvar)$LEGISLATIONR12)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~LEGISLATIONR12, ~stay_leave, dsgn.never19, svyvar)$LEGISLATIONR12)/sqrt(nrow(d19_Jo_Nv)))
# LEGISLATIONR13
y13 <- c(svyby(~LEGISLATIONR13, ~ stay_leave, dsgn.stay19, svymean)$LEGISLATIONR13,svyby(~LEGISLATIONR13, ~ stay_leave, dsgn.never19, svymean)$LEGISLATIONR13)
y13_error <- c(sqrt(svyby(~LEGISLATIONR13, ~stay_leave, dsgn.stay19, svyvar)$LEGISLATIONR13)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~LEGISLATIONR13, ~stay_leave, dsgn.never19, svyvar)$LEGISLATIONR13)/sqrt(nrow(d19_Jo_Nv)))
# LEGISLATIONR14
y14 <- c(svyby(~LEGISLATIONR14, ~ stay_leave, dsgn.stay19, svymean)$LEGISLATIONR14,svyby(~LEGISLATIONR14, ~ stay_leave, dsgn.never19, svymean)$LEGISLATIONR14)
y14_error <- c(sqrt(svyby(~LEGISLATIONR14, ~stay_leave, dsgn.stay19, svyvar)$LEGISLATIONR14)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~LEGISLATIONR14, ~stay_leave, dsgn.never19, svyvar)$LEGISLATIONR14)/sqrt(nrow(d19_Jo_Nv)))
# LEGISLATIONR15
y15 <- c(svyby(~LEGISLATIONR15, ~ stay_leave, dsgn.stay19, svymean)$LEGISLATIONR15,svyby(~LEGISLATIONR15, ~ stay_leave, dsgn.never19, svymean)$LEGISLATIONR15)
y15_error <- c(sqrt(svyby(~LEGISLATIONR15, ~stay_leave, dsgn.stay19, svyvar)$LEGISLATIONR15)/sqrt(nrow(d19_St_Lv)),sqrt(svyby(~LEGISLATIONR15, ~stay_leave, dsgn.never19, svyvar)$LEGISLATIONR15)/sqrt(nrow(d19_Jo_Nv)))
## combining calculations
df3 <- rbind.data.frame(cbind(x,y=y,y_error=y_error),cbind(x,y=y2,y_error=y2_error),cbind(x,y=y3,y_error=y3_error),
                        cbind(x,y=y4,y_error=y4_error),cbind(x,y=y5,y_error=y5_error),cbind(x,y=y6,y_error=y6_error),
                        cbind(x,y=y7,y_error=y7_error),cbind(x,y=y8,y_error=y8_error),cbind(x,y=y9,y_error=y9_error),
                        cbind(x,y=y10,y_error=y10_error),cbind(x,y=y11,y_error=y11_error),cbind(x,y=y12,y_error=y12_error),
                        cbind(x,y=y13,y_error=y13_error),cbind(x,y=y14,y_error=y14_error),cbind(x,y=y15,y_error=y15_error))
df3$y <- -1*as.numeric(as.character(df3$y))
df3$y_error <- -1*as.numeric(as.character(df3$y_error))
df3$measure <- c(rep("Whites",4),rep("Blacks",4),rep("Hispanics",4),rep("People who live in large cities",4),
                 rep("People in small towns/rural",4),rep("Poor people",4),
                 rep("Middle class people",4),rep("Wealthy people",4),rep("Legal immigrants",4),
                 rep("Illegal immigrants",4),rep("Young adults",4),rep("The elderly",4),
                 rep("Large corporations",4),rep("Colleges and universities",4),rep("Churches/places of worship",4))
colnames(df3) <- c("model","estimate","std.error","term")
# getting rank from other dataset and merging in for consistency
df2.rank <- df2[c("term","rank")]
df2.rank <- unique(df2.rank)
df3.m <- merge(df3, df2.rank, by = "term")
df3.m <- rbind.data.frame(df3.m[df3.m$rank==15,],df3.m[df3.m$rank==14,],df3.m[df3.m$rank==13,],df3.m[df3.m$rank==12,],
                          df3.m[df3.m$rank==11,],df3.m[df3.m$rank==10,],df3.m[df3.m$rank==9,],df3.m[df3.m$rank==8,],
                          df3.m[df3.m$rank==7,],df3.m[df3.m$rank==6,],df3.m[df3.m$rank==5,],df3.m[df3.m$rank==4,],
                          df3.m[df3.m$rank==3,],df3.m[df3.m$rank==2,],df3.m[df3.m$rank==1,])
## plotting data
# taxes plot
p1 <- dwplot(df2, dodge_size = 0.4, dot_args=list(aes(colour = model, shape = model), size = 3), whisker_args=list(aes(colour = model))) + 
  theme_minimal() + xlab("RECEIVE MORE in benefits than they pay in taxes --to-- RECEIVE LESS than they pay") + 
  theme(legend.title = element_blank(), legend.position = "top",legend.key.width = unit(.5,"cm"), 
        plot.title = element_text(size = 12, hjust=0.5), axis.title.x = element_text(size = 9), axis.text = element_text(size=10)) + 
  xlim(-.75,.75) + scale_colour_manual(breaks=c("Dem Loyalists","Dem Switchers","Rep Switchers","Rep Loyalists"), values= c("darkorange", alpha("darkorange", .45), alpha("purple3", .45), "purple3")) +
  scale_shape_manual(breaks=c("Dem Loyalists","Dem Switchers","Rep Switchers","Rep Loyalists"), values = c(20,20,20,20)) +
  labs(title = "Taxes")
# legislation plot
p2 <- dwplot(df3.m, dodge_size = 0.4, dot_args=list(aes(colour = model, shape = model), size = 3), whisker_args=list(aes(colour = model))) + 
  theme_minimal() + xlab("U.S. laws are TOO FAVORABLE --to-- Laws are NOT FAVORABLE ENOUGH") + 
  theme(legend.title = element_blank(), legend.position = "top",legend.key.width = unit(.5,"cm"), 
        plot.title = element_text(size = 12, hjust=0.5), axis.title.x = element_text(size = 9), axis.text = element_text(size=10)) + 
  xlim(-.75,.75) + scale_colour_manual(breaks=c("Dem Loyalists","Dem Switchers","Rep Switchers","Rep Loyalists"), values= c("darkorange", alpha("darkorange", .45), alpha("purple3", .45), "purple3")) +
  scale_shape_manual(breaks=c("Dem Loyalists","Dem Switchers","Rep Switchers","Rep Loyalists"), values = c(20,20,20,20)) +
  labs(title = "Legislation")
# combining plots
legend <- get_legend(p1)
p1 <- p1 + theme(legend.position="none")
p2 <- p2 + theme(legend.position="none")
g1 <- grid.arrange(legend, p1, p2, ncol=1,
                   layout_matrix = rbind(1,2,3),
                   widths = 4, heights = c(.25,2,2))
ggsave("online_appendix_A/plots/FigG_laws_taxes_items.png", g1, width = 7, height = 7)



